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The statistics of ACDM Halo Concentrations 
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ABSTRACT 

We use the Millennium Simulation (MS) to study the statistics of ACDM halo con- 
centrations at z = 0. Our results confirm that the average halo concentration declines 
monotonically with mass; a power-law fits well the concentration-mass relation for 
over 3 decades in mass, up to the most massive objects to form in a ACDM uni- 
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verse (~ 10 h M^). Thi s is in clear disagreement with the predictions of the model 



proposed by iBullock et al. for these rare objects, and agrees better with the original 
predictions of lNavarro, Frenk. fc Whitel The large volume surveyed, together with the 
unprecedented numerical resolution of the MS, allow us to estimate with confidence 
the distribution of concentrations and, consequently, the abundance of systems with 
unusual properties. About one in a hundred cluster haloes (M 200 >3 x 10 14 h~ x M Q ) 
have concentrations exceeding C200 = 7.5, a result that may be used to inter- 
pret the likelihood of unusually strong massive gravitational lenses, such as Abell 
1689, in the ACDM cosmogony. A similar fraction (1 in 100) of galaxy-sized haloes 
(M200 ~ 10 12 hr 1 M Q ) have C200 < 4.5, an important constraint on models that at- 
tempt to reconcile the rotation curves of low surface-brightness galaxies by appealing 
to haloes of unexpectedly low concentration. We find that halo concentrations are 
independent of spin once haloes manifestly out of equilibrium are removed from the 
sample. Compared to their relaxed brethren, the concentrations of out-of-equilibrium 
haloes tend to be lower and to have more scatter, while their spins tend to be higher. 
A number of previously noted trends within the halo population are induced primarily 
by these properties of unrelaxed systems. Finally, we compare the result of predicting 
halo concentrations using the mass assembly history of the main 'progenitor with esti- 
mates based on simple arguments based on the assembly time of all progenitors. The 
latter typically do as well or better than the former, suggesting that halo concentra- 
tion depends not only on the evolutionary path of a halo's main progenitor, but on 
how and when all of its constituents collapsed to form non-linear objects. 
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1 INTRODUCTION 

The advent of large cosmological N-body simulations has en- 
abled important progress in our understanding of the struc- 
ture of dark matter haloes. As a result, over the past few 
years a broad consensus has emerged about the mass func- 
tion of these collapsed structures, about their mass profiles 
and shapes, and about the presence of substructure within 
the vitalised region of a halo. The spherically-averaged halo 



mass profiles are of particular interest, not only because of 
their immediate applicability to a host of observational di- 
agnostics, such as gravitational lensing and disk galaxy ro- 
tation curves, but also because of the apparent simplicity of 
their structure. 

iNavarro. Frenk. fc Whitel \ 19951 . 1 19961 . Il997l . hereafter 
NFW) argued that the density profile of a dark matter halo 
may be approximated by a simple formula with two free 
parameters; 
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where p crlt — 3Hq/8itG is the critical density for closurcQ, 
5 C is a characteristic density contrast, and r 3 is a scale ra- 
dius. Remarkably, this formula seems to hold for essentially 
all haloes assembled hierarchically and close to virial equi- 
librium, regardless of mass and of the details of the cosmo- 
logical model. The cosmological information is encoded in 
correlations between the parameters of the NFW profile, so 
that observational constraints on such parameters may be 
translated directly into interesting constraints on cosmolog- 
ical parameters. 

As discussed by NFW, such correlations arise be- 
cause the characteristic density of a system appears to 
evoke the density of the universe at a suitably defined 
time of collapse. This result has been revisited and 
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which has led to the development of a number of semi- 
analytic and empirical procedures to explain and predict 
the structural parameters of cold dark matter (CDM) 
haloes as a function of mass, redshift, and cosmological 
parameters. The various approaches differ in detail and 
lead to significantly different predictions, especially when 
extrapolated to halo masses or to redshifts which were not 
well sampled by the numerical data on which they were 
based. 

NFW, for example, proposed that the characteristic 
density is set at the time when most of the mass of a halo 
is in non-linear, collapsed structures. iBullock et al.l l|200ll . 
B01), on the other hand, argued that a better fit to their N- 
body results is obtained by assuming that the scale radius 
of haloes of fixed virial mass is independent of redshift, lead- 
ing to a substantially different redshift evolution of the halo 
structural parameters th an envisioned by NFW. Th i s con- 
clusion was seconded by Ek e. Navarro, fc Steinmet j l|200ll . 
ENS) who proposed a modification of BOl's approach to take 
into account models with truncated power spectra, such as 
expected, for example, in a warm dark matter-dominated 
universe. 

The predictions of these models also differ significantly 
for extremely massive haloes, but these predictions have 
been notoriously difficult to validate since these systems 
are woefully under-represented in simulations that survey 
a small fraction of the Hubble volume. For example, NFW 
argued that the characteristic density of a halo is set by the 
mean density of the Universe at the time of collapse. Very 
massive systems have, by necessity, been assembled quite 
recently (indeed, they are assembling today), and therefore 
they should all have similar characteristic densities. BOl's 
model, on the other hand, defines collapse redshifts in a dif- 
ferent way, predicting much lower concentrations for very 
massive objects. 

Despite the (qualified) success of these models at re- 
producing the average mass and redshift dependence of halo 
structural parameters, they are in general unable to account 
for the siz able scatter about the mean relations. As first dis- 
cussed bv lJind l|2000h . a sizable spread in concentration (of 
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ber of models have attempted to reproduce this result us- 
ing semi-analytic models. Amongst the most successful are 
those that ascribe variations in concentration to disparities 
in the assembly histor y of haloes of given mass. For example, 
IWechsler et al.l (|2002l . W02) identify the scatter with varia- 
tions in the time when the rate of mass accretion onto the 
main progenitor pe aks. A similar proposal was advanced by 
I Zhao et al.l (|2003al fbl. Z01), who argued that the concentra- 
tion of a halo is effectively set during periods when the most 
massive progenitor is in a phase of fast mass accretion. 

Given the disparities between models, it is important to 
validate their predictions in a regime different from that used 
to calibrate their parameters. Indeed, with few exceptions, 
most of these studies have explored numerically a relatively 
narrow range of halo mass and redshift, favouring (because 
they are easier to simulate) haloes with masses of the or- 
der of the characteristic non- linear mass, M* , and redshifts 
close to the present day (z ~ 0). Testing these predictions 
on a representative sample of haloes of mass much greater 
or much lower than M,, or at very high redshift, requires 
either simulations of enormous dynamic range, or especially 
designed sets of simulations that probe various mass or red- 
shift intervals one at a time. 

One version of the latter approach was adopted by 
NFW, who simulated individually haloes spanning a large 
range in mass. The price paid is the relatively few haloes 
that can be studied using such simulation series, as well as 
the lingering possibility that the procedure used to select 
the few simulated h aloes may int r oduce some subtle bias 
or artifact. Recently, iMaccio et al. I (|2007l , M07), have com- 
bined several simulations of varying mass resolution and box 
size in order to try and extend earlier results to M < M* 
scales. This approach is not without pitfalls, however. For 
example, in order to resolve a statistically significant sam- 
ple of haloes with masses as low as 10 10 /i -1 Mq, M07 use a 
simulation box as small as 14.2 h -1 Mpc on a side, leading 
to concerns that the substantial large-scale power missing 
from such small periodic realisation may unduly influence 
the results. 

One way to overcome such shortcomings is to increase 
the dynamic range of the simulation, so as to encompass a 
volume large enough to be representative while at the same 
time having enough mass resolution to extend the analy- 
sis well below or well above M*. This is the approach we 
adopt in this paper, where we use the Millennium Simula- 
tion (MS) to address these issues. The enormous volume of 
the MS (500 3 /i~ 3 Mpc 3 ), combined with the vast number of 
particles (2160 3 ), make this simulation ideal to characterise, 
with minimal statistical uncertainty, the dependence of the 
structural parameters of ACDM haloes on mass, spin, for- 
mation time, and departures from equilibrium. We extend 
and check our MS results at low masses by using an addi- 
tional large simulation of a 100 ft -1 Mpc region with about 
10 times better mass resolution. 

For reasons discussed in detail below, numerical limita- 
tions impose a lower mass limit of about 10 12 /i -1 Mq in our 
analysis. Thus, our study does not extend to halo masses as 
low as those probed by M07 (whose smallest box is filled with 
particles of mass 1.4 x 10 7 h -1 Mq), but is aimed at extend- 
ing previous work to give reliable and statistically robust 
results for large and representative samples of haloes over 
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the full ranee from 10 12 to 10 15 /i _1 M . The large number 
of haloes in the MS also allows us to study in detail devia- 
tions from the mean trends and, in particular, the possible 
presence of systems with unusual properties, such as clus- 
ters with unusually high concentrations or galaxy haloes of 
unusually low density. In this paper we concentrate on the 
properties of haloes at z — 0. A second paper will extend 
these results to high redshift. 

The plan for this paper is as follows. We describe briefly 
the simulation and the halo identification technique in Sec- 
tion [2] After a brief introduction to the MS (Section l2.1[) we 
describe our halo identification (Section I2.2[l and selection 
(Section 12. 3|) techniques. We also describe the merger trees 
in Section ^. 4l and the NFW profile fitting procedure in Sec- 
tion 12.51 The dependence of the halo structural parameters 
on mass, spin, and formation time, as well as the perfor- 
mance of various semi-analytic models designed to predict 
halo concentrations, are discussed in Section|3] We conclude 
with a brief summary in Section [4] 



2 HALOES IN THE MILLENNIUM 
SIMULATION 

Our analysis is mainly based on haloes identified in the Mil- 
lennium Simulation (MS), (|Springel et al.l feOQSal ). and in 
this section we describe briefly our halo identification and 
cataloguing procedure. For completeness, we begin with a 
brief summary of the main characteristics of the MS, and 
then move on to a fairly detailed characterisation of the halo 
sample. Readers less interested in these technical details may 
wish to gloss over this section and skip to Section [3] where 
our main results are presented and discussed. 



2.1 The simulations 

The Millennium Simulation is a large N-body simulation of 
the concordance ACDM cosmogony. It follows N — 2160 3 
particles in a periodic box of Lbox = 500/i -1 Mpc on a 
side. The cosmological parameters were chosen to be consis- 
tent w ith a combined analy sis of the 2dFGRS (jColless et al.l 

year WMAP data 
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i Spergel et all [20031. They are 0, 
Q h = 0.045, h = 0.73, f2 A = 0.75, 
f2 denotes the present day contribution of each component 
to the matter-energy density of the Universe, expressed in 
units of the critical density for closure, p cr it; n is the spectral 
index of the primordial density fluctuations, and as is the 
linear rms mass fluctuations in 8/i _1 Mpc spheres at z = 0. 
Compared with the parame ter values now favou red by the 
three- year WMAP analysis (ISpergel et al. l |2006t), the main 
differences are that a modest tilt, n — 0.95 ± 0.02 and a 
lower erg = 0.74 ± 0.05 are favoured by the analysis of these 
latest data. 

With our choice of cosmological parameters, the par- 
ticle mass in the MS is 8.6 x 10 8 /i _1 Mq. Particle pairwise 
interactions are softened on scales smaller than (Plummer- 
equivalent) e = 5/i~ 1 kpc. Since galaxy-sized haloes (M ~ 
10 12 hT 1 Mq) in the MS are represented with only about 
1000 particles, we have verified that our results are insensi- 
tive to numerical resolution by comparing them with a sec- 
ond simulation of a smaller volume, Lbox = 100 fo _1 Mpc, but 



of 9 times higher mass resolution. This simulation adopted 
the same cosmological model as the MS, and evolved N = 
900 3 particles of mass 9.5 x 10 7 /i _1 Mq, softened on scales 
smaller than e = 2.4/i _1 kpc. 

Both simulations were performed with a special ver- 
sion of the GADGET-2 code (|Springelll2005bl l that was spe- 
cially designed for massively parallel computation and for 
low memory consumption, a prerequisite for a simulation of 
the size and computational cost of the MS. 



2.2 Halo identification 

The simulation co de produced on the fly a friends of friends 
jDavis et al.lll985l ) (FOF) group catalogue with link param- 
eter, b = 0.2, and at least 20 particles per group. At z = 0, 
this procedure identifies 1.77 x 10 6 groups in the MS. We 
have also used SUBFIND , the s ubhalo finder algorithm de- 
scribed in ISpringel et ail l|200ll l. in order to clean up the 
group catalogue of loosely-bound FOF structures, and to 
analyse the substructure within each halo (see Sec. 12. 3[) . 

Like FOF, SUBFIND keeps only substructures contain- 
ing 20 or more particles. In this way, each FOF halo is 
decomposed into a background halo (or the most massive 
"substructure") and zero or more embedded substructures. 
In the MS, SUBFIND finds at z = a total of 1.82 x 10 7 
substructures, with the largest FOF group containing 2328 
of them. 



2.2.1 Halo centring 

Since much of our analysis deals with radial profiles, it is im- 
portant to define carefully the centre of a halo. We choose 
the position, r c , of the particle with minimum gravitational 
potential in the most massive substructure (the potential 
centre). Although this seems like a sensible choice, it is im- 
portant to check that other plausible options do not lead to 
large differences in the location of the halo centre. 

We have therefore compared the potential centre with 
the r esult of the "shrinking sphere" algorithm l|Power et al.l 
2003), which is intended to converge towards the density 
maximum of the most massive substructure, independent 
of the SUBFIND algorithm. It starts by enclosing all FOF 
particles within a sphere and computes iteratively their cen- 
tre of mass, shrinking the radius of the sphere by r* = 
ro(l — 0.025) ! , and rejecting particles outside the sphere. 
The iteration stops when the shrinking sphere contains 1% 
of the initial number of particles. 

We carried out this comparison in a sub-volume of the 
MS containing 2000 haloes with Afof > 450. For 93% of 
these haloes the methods agree in the centre position, with 
a difference smaller than the gravitational softening, e. How- 
ever, we note that the result can depend on the geometry of 
the FOF group. When the FOF halo is double (or multiple) 
and its centre of mass is far from the centre of the most 
massive substructure, the shrinking sphere may converge to 
another slightly less massive substructure. We conclude that 
the potential centre is a more robust determination of the 
halo centre, but that the discrepancy between the two meth- 
ods could be used to flag problematic haloes whose mass 
distributions deviate significantly from spherical symmetry. 
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2.2.2 Halo boundary 



Using the potential centre, we define the limiting radius rn m 
of a halo by the radius that contains a specified density 
contrast ~p(r) = Ap cr it. This defines implicitly an associated 
mass for the halo through 
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We note that this includes all the particles inside this spher- 
ical volume, and not only the particles grouped by the FOF 
or the SUBFIND algorithms. 

The choice of A varies in the literature, with some 
authors using a fixed value, such as NFW, who adopted 
A = 200, and others, such as B01, who choose a value moti- 
vated by the spherical collapse model, where A ~ f 78 fi^ 45 
(for a fiat universe), which gives A = 95.4 at z = 
for our adopted ACDM p arameters (e.g. lLahav et al.lll99ll ; 
lEke. Cole, fc Frenklll996l l . The drawback of the latter choice 
is its dependence on redshift and cosmological parameters. 
We keep track of both definitions in our halo catalogue, but 
will quote mainly values adopting A = 200. When necessary, 
we shall specify the choice by a subscript; e.g., M200 and r2oo 
are the mass and radius of a halo adopting A = 200; M v i r 
and r v j r correspond to adopting A = 95.4. Unless otherwise 
specified, quantities listed without subscript throughout the 
paper assume A = 200. 



2.3 Halo selection 

Dark matter haloes are dynamic structures, constantly ac- 
creting material and often substantially out of virial equi- 
librium. In these circumstances, haloes evolve quickly, so 
that the parameters used to specify their properties change 
rapidly and are thus ill-defined. Furthermore, in the case of 
an ongoing major merger, even the definition of the halo 
centre becomes ambiguous, so that the characterisation of a 
system by spherically-averaged profiles is of little use. As we 
shall see below, departures from equilibrium not only add 
to the scatter in the correlations that we seek to establish, 
but can also bias the resulting trends, unless care is taken 
to identify and correct for the effect of these transient struc- 
tures. 



2.3.1 Relaxed and unrelaxed haloes 

The equilibrium state of each halo is assessed by means of 
three objective criteria: 



i) Substructure mass fraction: We compute the mass 
fraction in resolved substructures whose centres lie inside 
r v iT- /sub = X/i5n b Msub,i/M V if Note that in this definition 



/sub does not include the most massive substructure as this 
is simply the bound component of the main halo. 

ii) Centre of mass displacement: We define s, the nor- 
malised offset between the centre of mass of the halo (com- 
puted using all particles within r v j r ) and the p otential cen- 
tre, as s = |r c — r cm |/r v i r l| Thomas et al.|[200lh . 

iii) Virial ratio: We compute 2T/\U\, where T is the total 
kinetic energy of the halo particles within r v i r and U their 
gravitational self potential energy. To estimate U, we use a 
random sample of 1000 particles when Ni > 1000. We obtain 
physical velocities with respect to the potential centre by 
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Figure 1. Images (top) and corresponding spherically averaged 
density profiles (bottom) in four haloes of similar mass. The halo 
shown in the lower right panel of each set satisfies all our se- 
lection criteria and is, therefore, close to dynamical equilibrium. 
Note that the NFW profile (solid line) provides an excellent fit to 
this halo. The halo mass, concentration and the values of the three 
quantities, / su b, s and 2T/\U\ used in the selection are given in 
the legend. The NFW fitting procedure, which here is performed 
only over the indicated range r m i n < r < r v i r , is described in 
Section 12.51 The remaining three haloes are excluded from our 
relaxed sample as they fail at least one of the selection criteria. 
The halo on the upper left has a large amount of substructure, 
/sub > 0.1. The one in the upper right panel is undergoing a ma- 
jor merger. Note that the merging partner does not contribute to 
/sub) since its centre lies outside the virial radius, but some of its 
associated material displaces the centre of mass of the system, re- 
sulting in s > 0.07. The halo in the lower left panel satisfies these 
two criteria, but has 2T/\U\ > 1.35. The corresponding panels 
in the lower plot show that these unrelaxed haloes have density 
profiles that are clearly not well described by NFW profiles. 
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adding the Hubble flow to the peculiar velocities and then 
we compute the halo kinetic energy after subtracting from 
the velocities the motion of the halo centre of mass. 

Fig. [1] shows images and spherically averaged density 
profiles for a set of three haloes of similar mass that are re- 
jected by just one of each of the above criteria. Systems such 
as the one shown in the top-left panel (large / au b) are clearly 
not relaxed and it is not surprising that the large number of 
substructures affect halo prop erties such as concentration, 
angular momentum and shape (|Gao et al. 2004; Sha w et al.l 
120051 ^1 . Note that, despite the large value of / su b, the centre 
of mass displacement is quite small, since the spatial distri- 
bution of substructures happens to be fairly symmetric. 

Conversely, the halo shown in the upper-right panel of 
Fig.[l]has small f su b, but is rejected by our cut on the centre 
of mass offset s. This is clearly an ongoing merger where, be- 
cause the merging partner of the main halo lies just outside 
the virial radius, it contributes little to /sub- Thus the s cri- 
terion is complementary to / su b- This is important, because 
the ability of SUBFIND to detect self-bound substructures 
is heavily resolution-dependent: / su b is likely to be substan- 
tially underestimated in low-mass haloes resolved with few 
particles. 

These two criteria, however, make no use of kinematic 
information and so our third criterion, based on 2T/\U\, 
is a useful supplement able to reject haloes that, despite 
passing the other two conditions, are far from dynamical 
equilibrium. This is especially the case for ongoing mergers 
and artificially linked haloes. An example of a halo rejected 
by just this criterion is shown in the lower-left panel of Fig. [T] 

Finally, the lower- right panel of Fig.[T]shows an example 
of a halo that makes it into our relaxed sample. These relaxed 
haloes generally have smooth density profiles and are well 
fitted by NFW profiles. 



2.3.2 N-dependence of selection criteria 

Fig. [5] shows the various correlations between the criteria 
used to identify "relaxed" haloes (/ su b, s and 2T /\U\) and 
the number of particles within the virial radius. Note that 
the equilibrium measures we use do not have bimodal distri- 
butions, but, instead, are roughly continuous, with extended 
tails. This reflects the fact that haloes assemble hierarchi- 
cally: all haloes are in the process of accreting some material 
and are, therefore, to some extent unrelaxed. (Note for ex- 
ample, that the median 2T/\U\ is slig htly greater than unit y, 
as reported in earlier work (see, e.g. ICole fc Lacev|[T996t) ). 
Thus, our selection criteria just provide a simple, but some- 
what arbitrary, way of trimming off the tail of worst offend- 
ers, typically ongoing mergers. 

The shaded regions in Fig. [2] show the areas of param- 
eter space rejected by our criteria to select relaxed haloes: 
/ sub < 0.1, s < 0.07 and 2T/\U\ < 1.35. The top-right panel 
shows the expected strong correlation between our measures 
of asymmetry and lumpiness, as well as how the tail of high 
s and high / su b values is removed by the selection. 

Of all the selection criteria, the one most sensitive to 
the number of particles in the halo is / au b. This is clearly 
seen in the top-left panel of Fig. [2] where the stepped line 
shows, as a function of iV v j r , the fraction of haloes with 
no resolved substructures, i.e., / su b = 0. This rises quickly 
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Figure 2. The various criteria used to define our relaxed sample 
of haloes (see Section |2. 31 for definitions), shown as a function of 
the number of particles within the virial radius. Shaded regions 
indicate the location of "unrelaxed" haloes. The dots sample uni- 
formly all MS haloes with AT vlr > 300 in each plot. Top left: The 
fraction of halo mass in resolved substructures, / su b vs. the num- 
ber of particles in the halo N. The dashed line shows the detection 
limit of SUBFIND, Af min = 20/N. The solid stepped line shows 
the fraction of haloes with no resolved substructures in each mass 
bin (i.e., those with / su b = 0). Top right: The centre of mass offset 
s vs the substructure fraction / su b- The criteria / su b < 0.1 and 
s < 0.07 reject the tail of haloes with high / su b and s. Bottom: 
The criterion 2T/\U\ < 1.35 is intended to reject haloes far from 
dynamical equilibrium. 



with decreasing JV v ir, so that more than 80% of haloes with 
300 < iVvir < 1000 particles have no discernible substruc- 
tures. By contrast, essentially all haloes with more than 
10, 000 particles have at least one massive subhalo within 
the virial radius. 

A significant fraction of haloes with fewer than 1000 
particles also have very large values of the virial ratio 
2T/\U\, a s shown in t he bo ttom-left panel of Fig. [2] As dis- 
cussed bv lBett et ail (|2007[ ). these are loosely bound objects 
connected by a tenuous bridge of particles, or objects lying 
in the periphery of much more massive systems that owe 
their large kinetic energy to contamination with fast-moving 
members of the other system. 

These results suggest that a minimum number of par- 
ticles may also be required in order to obtain robust results 
independent of numerical artifact. We shall see below that 
about ~ 1000 particles or more are needed in order to avoid 
biases when deriving structural parameters from fits to the 
halo mass profiles. 

The fraction of haloes rejected by these cuts varies 
slowly with mass, rising from 20% at M v i r = 10 12 M© to 
50% at M vir = 10 15 h~ 1 M e (see numbers given in Fig.|6]and 
Table [TJ. The criterion s < 0.07 rejects the most haloes; a 
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Figure 3. Density profiles, r 2 p(r), and least squares NFW fits for 
four relaxed haloes. The fits are performed over the radial range 
0.05 < r/r v j r < 1, shown by the solid circles, and extend slightly 
beyond r2oo- The vertical line marks the position of the mass 
profile convergence radius derived by Power et al (2003). This is 
always smaller than the minimum fit radius, 0.05 r v i r . 



smaller but significant fraction of haloes that pass that cut 
are removed by the / su b < 0.1 criterion; while only a few 
additional haloes are rejected by the 2T/\U\ < 1.35 test. 



2.4 Merger trees and formation times 

We use the merger trees for the MS d e scribe d in detail by 
iHarker et all (|2006l ) and iBower et al] d2006l ). These differ 
slightly from those used bv lSpringell ( 2005bh . but the differ- 
ences are only significant for haloes undergoing major merg- 
ers, which would not pass the stringent selection criteria 
listed above. At each of the approximately 60 output red- 
shifts of the simulation, the merger tree provides us with a 
list of all the haloes that will subsequently merge to become 
part of the final halo. We exploit this information below in 
order to investigate the dependence of halo concentration on 
formation history (see Section f3 . 4 p . 



2.5 Profile fitting 

For each halo identified using the procedure outlined in Sec- 
tion 12.21 we have computed a spherically-averaged density 
profile by binning the halo mass in equally spaced bins in 
log 10 (r), between the virial radius and log 10 (r/r v i r ) = —2.5. 
After extensive testing, we concluded that using 32 bins, 
corresponding to Alog 10 (r) = —0.078, is enough to produce 
robust and unbiased results. 

The density profiles of four relaxed haloes are shown in 
Fig. El together with fits using the NFW profile (eq.HJ. This 
profile has two free parameters, 8 C and r S) which are both 
adjusted to minimise the rms deviation, (Tat, between the 
binned p(r) and the NFW profile, 

1 JVbinS 2 

<7fit = jj — r V [log 10 pi - log 10 p N Fw(<5 c ; r s )] . (3) 

-^bins r 

i=l 

Note that eq. [3] assigns equal weight to each bin. We have 



explicitly checked that none of our results varies significantly 
if we adopt other plausible choices, s uch as Poiss on weighting 
each bin (for further discussion, see I Jind (|2000l )). 

Once the parameters 8 C and r s for each halo are re- 
trieved from the fitting procedure, they may be expressed in 
a variety of forms. In order to be consistent with the original 
NFW work, we choose to express the results in terms of a 
mass, M200, and a concentration, c = C200 = r 2oo/r s , that 
result from adopting A = 200 in the definition of the virial 
properties of a halo (eq.[2]). 

The radial range adopted for the fitting procedure is im- 
portant, especially since haloes of different mass are resolved 
to varying degree in a single cosmological simulation. After 
experimenting at length, and especially after comparing the 
fit parameters obtained in the overlapping mass range of the 
two simulations ( Section [2TT}, we settled on carrying out the 
fits over a uniform radial range (in virial units) . This ensures 
that all haloes, regardless of mass, are treated equally, min- 
imising the possibility of introducing subtle mass-dependent 
biases in the analysis. 

Fig |3] shows the mass-concentration dependence ob- 
tained for the MS (symbols) and the higher mass resolu- 
tion simulation (lines) for four different choices of the radial 
range. The symbols and lines extend down to haloes with 
~ 1000 particles within the virial radius in each case and 
indicate the median (solid circles and lines) as well as the 
20 and 80 percentiles of the distribution. To guide the com- 
parison, the dotted line shows a power law, c oc M^j/ 10 , and 
is the same in all panels. 

Note that, provided the minimum radius is > 0.05 r v i r , 
the fit parameters appear to depend only weakly on the 
radial range, but that the distribution of concentrations is 
narrowest for 0.05 < r/r v i r < 1. For this choice (which we 
adopt hereafter in our analysis) there is also good agreement 
between the two simulations on mass scales probed in both 
with at least 1000 particles. 

However, as shown in Fig. [5] the mean quality of the fits, 
as measured by Cfi t , deteriorates noticeably for haloes with 
less than 10, 000 particles as a result of the relatively small 
number of particles per bin. Although this does not seem 
to introduce a bias in the mass-concentration relation, we 
have decided to retain, conservatively, only haloes with N > 
10, 000 for our analysis, and to combine the two simulations 
in order to probe the halo mass range 10 12 < M/h~ 1 Mo < 
10 15 . 

Finally, we considered a couple of alternative methods 
for characterising the halo concentration that do not as- 
sume an NFW density profile, such as th e ratio of maxi - 
mum to virial circular velocities, Vm^/V v - lr l|Gao et al.ll2004l ) 
or the ratio of masses enclosing diff erent overdensities, such 
as M v i r /Miooo (|Thomas et alj|200ll ). These measures of the 
concentration, presumably because they use information on 
just two particular points in the profile, lead to substantially 
larger scatter in the correlations that we examine here, so 
we do not consider them further in this paper. 



3 RESULTS AND DISCUSSION 
3.1 Concentration vs mass 

Fig. [6] shows concentration as a function of halo mass for 
all haloes selected following the procedure described above. 
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Figure 4. The median, 20 and 80 percentiles of the concentration, 
C20O1 as a function of halo mass, M2oo- The symbols extending to 
high masses show the results from the MS while the overlapping 
set of solid and dotted lines extending to lower masses show the 
results from a simulation with 9x higher mass resolution. Each 
panel corresponds to a different radial range adopted for the fits. 
Data for each simulation are shown for haloes with N > 1000 
particles, corresponding to ~ lO 12 /^ 1 Mq in the MS. The dotted 
line shows a power law, c <x M~ 1 ' 10 , and is the same in all panels. 



The upper panel is for our relaxed halo sample, while the 
lower panel show results for the complete sample, including 
systems that do not meet our equilibrium criteria. 

In both samples, the correlation between mass and con- 
centration is well defined, but rather weak. A power-law fits 
the median concentration as a function of mass fairly well; 
we find: 

c 200 = 5.26 (A/ 2 oo/10 14 /- :! " 1 " J " 



Mg 



(4) 



for relaxed haloes, and 



C200 = 4.67 (A/2oo /10 14 ft" 1 M ) ' (5) 

for the complete halo sample. 

Our power-law fit is in good agreement with the results 
of M07, who find C200 ~ 5.6 (Maoo/lO 14 ^ 1 Mq)" ' 098 for 
the average concentration of their sample of relaxed haloes. 
These authors also report that concentrations are system- 
atically lower when considering the full sample of haloes. 
The small difference between our results and M07's may 
be due to variations between mean and median, as well as 
on the different criteria used to construct the relaxed halo 
sample. Nevertheless, the agreement in the exponent of the 
power-law is remarkable, especially considering that these 
authors explore a mass range different from ours, namely 
2 x 10 9 < A/200 /ft - 1 A/ < 2 x 10 13 . Combining these re- 
sults with ours, we conclude that a single power law fits 
the concentration-mass dependence for about six decades in 
mass. 

Over the mass range covered by our simulations the 
concentration-mass dependence is in reasonable agreement 
with the predictions of NFW and of ENS, as shown, re- 
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Figure 5. The dependence of the rms residual deviation, crflt, 
about the best-fitting NFW density profile on the number of par- 
ticles per halo and on halo concentration. The boxes show the 
medians and the 25% and 75% centiles of the distribution, while 
the whiskers show the 5% and 95% tails. The numbers along 
the top of each panel indicate the number of haloes within each 
bin. Top panels include all MS haloes with AT v it > 450 as no 
other selection criteria have been applied. The upturn in (jfj t for 
low-concentration relaxed haloes is due to the inclusion of haloes 
with less than 10, 000 particles. This upturn disappears once the 
N > 10, 000 criterion is imposed. 



spectively, by the dotted and dashed lines in Fig. . The 
agreement, however, is not perfect, and both models appear 
to underestimate somewhat the median concentration at the 
low mass end. 

At the high-mass end, where there is a hint that concen- 
trations are approaching a constant value, the NFW model 
does slightly better than ENS. This is because a constant 
concentration for very massive objects is implicit in the 
NFW model, but not in ENS nor in the model of B01, which 
is shown by a dot-dashed line in Fig. [U Both ENS and B01 
predict a strong decline in concentration at the very high 
mass end. For the parameters favoured by B01, the disagree- 
ment for A/ > 10 13 5 ft~ 1 A/Q is dramatic, and cautions, as al- 
ready pointed out bv lZhao et al.l (|2003bl ). against using this 
model for predicting the concentrations of massive haloes. 

Finally, we note that M07 argue that the B01 model re- 
produces their results better than ENS for haloes of mass a 
few times 10 9 ft _1 Mq. However, the differences between the 
two models only become appreciable below ~ lO lo ft _1 M , 
which corresponds to only about 700 particles in their 
highest-resolution simulation. Given (i) the large scatter in 

2 These predictions use the original parameters in those papers 
and have not been adjusted further, except for adopting the 
power-spectrum "shape" parameter, T = 0.15, as the best match 
to the power spectrum adopted for the MS. 
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Figure 6. Mass-concentration relation for relaxed haloes (top 
panel) and for all haloes (bottom panel). The boxes represent the 
25% and 75% centiles of the distribution, while the whiskers show 
the 5% and 95% tails. The numbers on the top of each panel indi- 
cate the number of haloes in each mass bin. The median concen- 
tration as a function of mass is shown by the solid line and is well 
fit by the linear relations given in the lege nd of each panel. The 
dot-dashed line shows the prediction of the iBullock et al.l (|200lh 
model (using cvir2.f available from the authors); the dashed and 
dotted lines corresponds to the Eke, Navarro, & Steinmctz (2001) 
and NFW models, respectively, with T = 0.15, as this approxi- 
mates best the input power spectrum of the MS. See text for 
further details. 



the correlation; (ii) our finding that at least 1000 particles 
are needed to produce unbiased results, and (iii) the worries 
expressed in Section [1] about the rather small box used by 
M07 to resolve low-mass haloes, we conclude that it would 
be premature to judge conclusively on the superiority of one 
model over the other at the low mass end. More definitive 
tests are likely to come either from much higher resolution 
simulations or from extending this analysis to higher red- 
shifts, where the two models predict different behaviour. 
This is a topic that we intend to address in a forthcoming 
paper (Gao et al. 2007, in preparation). 



The large volume covered by the MS means that we have 
a large number of haloes in our sample, even in the most 
massive bins of Fig. [6j corresponding to rare, very massive 
objects with masses approaching 10 15 /i _1 Mq. We use this 
to check the common assumption that, at given halo mass, 
concent rations are distributed following a lognormal distri- 
bution (|jindl2000l ). 

We examine this in Fig. [7] where we show the distribu- 
tion of concentrations in two mass bins. In each panel the 
top histogram corresponds to all haloes in our sample, and 
the lower histogram only to those deemed relaxed by the 
criteria listed in Section \2. 31 

Although the overall distribution is not approximated 
well by a lognormal function, that of relaxed haloes is. The 
thick lines in Fig. [7] show fits of the form 



P(lo gl 



exp 



4( 



logio c - ( lo Sio c > 



(6) 



which is clearly a good approximation for the concentration 
distribution of relaxed haloes. 

Intriguingly, the concentration distribution of "unre- 
laxed" haloes may also be approximated by a lognormal 
distribution, albeit of smaller mean and larger dispersion, 
as shown by the dotted curves in Fig. [7] The sum of the 
two lognormals (shown with dashed lines) is indeed a very 
good approximation to the overall distribution. We list the 
median concentration and dispersion of these distributions, 
as well as the fraction of "unrelaxed" haloes, as a function 
of mass in Table Q] 

Upon inspection, a weak but systematic trend is appar- 
ent in Table [T] the dispersion in concentration appears to 
decrease monotonically as a function of mass (see also the 
bottom panel of Fig[8)l . This suggests that massive haloes are 
somehow a more homogeneous population than their lower 
mass counterparts, and may reflect the fact that massive 
haloes are rare objects that have collapsed recently, whereas 
less massive systems have a much wider distribution of as- 
sembly redshifts. We shall return to this topic below. 

One may apply the distributions shown in Fig [7] to esti- 
mate the abundance of relatively rare objects. These might 
be, for example, massive cluster haloes of unusually high 
concentration that stand out in X-ray or lensing surveys, or 
perhaps galaxy-sized haloes of unusually low concentration 
which may be notable because of peculiarities in the rotation 
curves of the galaxies they may host. 

We show this in the top and middle panels of Fig. [8] 
where the "plus" signs denote the fraction of haloes in each 
mass bin that have c < 4.5. The diamond symbols, on the 
other hand, indicate the fraction of systems with concentra- 
tions exceeding 7.5. For example, we find that slightly more 
than 1% of relaxed haloes with M200 ~ 3 x 10 14 /i _1 Mq have 
c > 7.5, or about 4 objects in our 500 3 fe~ 3 Mpc 3 volume. So 
clusters such as Abell 1689, for which iLimousin et all (|2006h 
estimate C200 = 7.6 should not be very abundant in the local 
Universe. Finding, through lensing or X-ray studies, many 
clusters as massive and as concentrated as Abell 1689 would 
yield strong constraints on the viability of ACDM as a cos- 
mological model. 
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Figure 7. The distribution of concentrations for haloes in two 
different mass bins: 11.75 < log 10 M20o/h~ 1 Mq < 12.25 (top 
panel) and 14.25 < log 10 M 2 oo//i _1 Mq < 14.75 (bottom panel). 
The top histogram in each panel shows the distribution for all 
haloes, the lower histogram that corresponding to the relaxed halo 
sample. The smooth curves are lognormal fits. Note that the over- 
all distribution may be very well approximated by the sum of two 
lognormal functions with parameters listed in Table [T] 



3.2 Concentration vs spin 

Although models such as those proposed by NFW, B01 or 
ENS may more or less account for the weak correlation be- 
tween mass and concentration, they say little about the ori- 
gin of the sizable scatter about the mean trend. As shown 
in Fig. [7] haloes of given mass have concentrations that may 
differ by up to a factor of three or more. What is the origin 
of this large scatter? 

One possibility is that concentration may be related 
to the angular momentum of the halo and that the large 
scatter in the sp in parameter at given mass (see, e.g., 
iBett et al] l|2007l ) and references therein) is somehow re- 
lated to the spread in concentration. This was explored 
by NFW and B01, who concluded that there was no ob- 
vious correla t ion b etween spin and concentration; however, 
Baili n~et al.l ((2005) have recently argued that such cor- 
relation indeed exists, and speculate that this may ex- 
plain why low surface-brigh tness galaxies (LSBs) appear 
to inhabit low-density haloes l)de Blok. Bosma. fc McGauehl 
l2003l ; lMcGaugh, Barker, fc de Blokll2003l ). 

We revisit this issue in Fig. |9j where we show the spin 
parameter A = J\E\ 1/2 /GM h/2 (where J is the angular 
momentum with respect the potential centre and E is the 
binding energy of the halo) versus concentration for our 
full halo sample (top panel) and for relaxed haloes (bottom 
panel). Although the full halo sample shows a well-defined 
tail of low-concentration, fast-spinning haloes reminiscent of 
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Figure 8. The fraction of relaxed haloes (top panel) and all haloes 
(middle panel) with c < 7.5 (plus sign symbols) and c > 4.5 
(diamond-like symbols). The points connected by solid lines are 
the direct measurements, in mass bins. Error bars denote 1/%/jV 
uncertainties. The lower panel shows the measured rms scatter in 
log 10 c for all and relaxed haloes, respectively. 

iBailin et all 's claim, this essentially disappears when relaxed 
haloes are considered. 

As discussed above, the low concentrations of unrelaxed 
haloes measured by profile-fitting algorithms are a transient 
result of the rapidly evolving mass distribution that accom- 
panies an accretion event, and do not indicate a halo with 
lower-than-average density. 

However, this observation does not explain why, as is 
apparent from Fig. [9] unrelaxed systems also tend to have, 
on average, higher spi ns than their relaxed coun terparts. 
As argued recently by iD'Onghia fc Navarro! (|2007l ). this is 
likely due to the redistribution of mass and angular momen- 
tum that occurs during mergers, combined with the arbi- 
trary virial boundaries used to define a halo and to compute 
its spin. Accretion events bring mass and angular momen- 
tum into a system, but redistribute it in such a way that 
high- angular momentum material ends up preferentially on 
weakly bound orbits that may take it beyond the nominal 
virial radius, thereby reducing A. This effect seems to be par- 
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Figure 9. The dependence of the spin parameter on concentra- 
tion for the full halo sample (top panel) and for relaxed haloes 
(bottom panel). The iso-density contours enclose 65, 95 and 99% 
of the distribution while the individual points show the distribu- 
tion for the remaining 1% of the distribution. 



ticularly important during mergers, and leads to an overall 
reduct i on of A as the merger progresses (see also iGardneJ 
l|200ll ): IVitvitska et all l|2002h ). 

We conclude, in agr eement with M07, that the bulk 
of the effect reported by iBailin et al.l (|2005t ) is due to the 
inclusion of out-of-equilibrium haloes in their sample. A very 
weak correlation between c and A may be visible in the tilt 
of the contours in the relaxed halo panel of Fig. [9] but it 
seems too weak to have strong observational implications. 



3.3 Concentration vs formation time 

The mass-concentration dependence was originally inter- 
preted by NFW as the result of the dependence of halo for- 
mation time on mass: the halo characteristic density just 
reflects the mean density of the Universe at the time of col- 
lapse. Therefore, low-mass haloes are more concentrated be- 
cause they collapsed earlier, when the Universe was denser. 



NFW also showed that the bulk of the scatter in the con- 
centration at given mass is due to variations in formation 
redshift. In this section we revisit the NFW analysis using 
our simulations. 

The first thing to note is that there is no unique 
way of defining when a particular halo formed, and 
so different definitions have been adopted (lLacev fc Cold 
19931; iNavarro, Frenk. fc White! [l997l ; IWechsler et al.l 120021 : 



Zhao et al Il2003al ). The simplest and most widely used def- 



inition is to set the formation redshift of a halo as the time 
when the most massive progenitor first exceeds half of the 
final halo mass, M(zf) = M(z = 0)/2. We use this as our 
default definition, but we will also consider other variants 
below. 

Fig. [10] shows the dependence of the formation red- 
shift (expressed as 1 + Zf) on halo mass. As expected, the 
median formation redshift clearly declines with increasing 
mass, and the solid line in the figure shows a least-squares 
fit to the median concentration. Dots show a random selec- 
tion of haloes coloured according to their concentration. The 
gradual change in colour from top to bottom indicates that 
concentration and formation time are closely related. 

We investigate this further in Fig[TT] where we plot the 
offset in concentration from the mean mass-concentration 
relation (Fig.|S} versus the formation redshift offset from the 
mean mass- formation redshift relation (Fig. llOp for two mass 
bins. The strong correlation between residuals indicates that 
the scatter in concentration at fixed mass is mostly due to 
variations in formation time. 

The fraction of the scatter in concentration shown in 
Fig. [7] <t igM , accounted for by z/ variations may be esti- 
mated by comparing it with the rms scatter, o"i gz f, about 
the 1:1 line in Fig. 1111 The fractional reduction in the vari- 
ance, lofgM — °f g zf I /°fgM i s 35% for the lowest mass range, 



10 11 ' 3 < M 2 oo//i _1 M < 10 12 - 25 , and 12% for the highest 



-)12.25 

mass range, 10 i42E> < M 20 o/h~ 1 M Q < 10 i4r \ 

This implies that the scatter in concentration is not 
fully explained by differences in formation time alone, and 
that additional effects are at work. One possibility is that 
the addition al scatter is an environ mental effect, as recently 
proposed bv lWechsler et all ([2006 ). However, these authors 
find that the effect is restricted to low-mass haloes, whereas 
our results show additional scatter for high mass haloes as 
well. 

Another possibility is that our formation time defini- 
tion should be revised. After all, the fraction of halo mass 
enclosed within the scale radius, r s (which defines c), is much 
less than one-half of the halo mass for the typical concentra- 
tions shown in Fig. [6] Or, finally, it might be that the scatter 
is driven by some aspect of the merger history that is not 
fully captured by our default definition of Zf , which depends 
on a single component of the halo (its most massive pro- 
genitor) rather than on the full spectrum of fragments that 
coalesce to form the final halo. 

We explore this in Fig. 1121 where we compare the NFW 
characteristic density, S c , with four different definitions of 
the formation redshift. We use S c in this figure because, 
besides being equivalent to the concentration, c, it is eas- 
ier to interpret. Indeed, if the characteristic density really 
tracks the mean density of the universe at the time of for- 
mation, one would expect it to follow the "natural" scaling, 

S c OC (1 +Z{) 3 . 
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Figure 10. Halo formation time as a function of mass for re- 
laxed haloes. Formation times are defined as the time when the 
most massive progenitor reaches half of the final halo mass. The 
numbers along the top of the main panel indicate the number of 
haloes in each bin. The straight line is a least square fit to the me- 
dian concentration. The colour of the plotted points encodes the 
concentration as indicated by the colour bar. The gradual change 
in colour, from green in the upper region of the plot to red at 
the bottom, shows qualitatively that concentration depends sen- 
sitively on the formation time. See also Fig. 1111 



Three of the Zf definitions explored in Fig.[T2]are varia- 
tions of the one adopted above, and indicate the time when 
the most massive progenitor first exceeds 10, 25, and 50% 
of the final halo mass, Mo = M(z = 0). These correspond 
to the panels labelled M o /10, Mo/4, and M /2 in Fig. [H 
respectively. The last definition, on the other hand, follows 
the prescription of NFW, and identifies the time when the 
combined mass of all M > Mo/10 progenitors exceeds Mo/2. 

As is clear from Fig. 1121 the tightest relation around the 
"natural scaling" (shown as solid lines) corresponds to the 
NFW definition. This suggests that the full mass spectrum 
of clumps that assembles into the halo plays an important 
role in the final halo concentration and not just the most 
massive progenitor. 

We also note that in each panel of Fig. [12] it is the first 
(lowest redshift) bin that has the largest scatter and that is 
furthest form the "natural scaling" relation. This bin cor- 
responds to haloes that have been assembled very recently, 
and may contain haloes that, despite our selection criteria, 
are still unrelaxed and for which the structural parameters 
are ill defined. 



3.4 Concentration predictions 

It is clear from the above discussion that accurate predic- 
tions of the concentration require some knowledge of the 
halo's assembly history. This cautions against the common 
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Figure 11. Concentration offsets from the mean mass- 
concentration relation (Fig. |6j versus formation-time offsets from 
the mean mass-formation time relation (Fig. llOU for two mass 
bins, as labelled in each panel. The strong correlation between 
residuals implies that much of the scatter in the M-c relation 
is due to variations in the formation time of haloes of given 
mass. Boxes indicate the median and (25,75) percentiles, whiskers 
stretch to the (5,95) percentiles. The number of haloes in each 
Alog 1 g(l + Zf) bin is given in the legend. 



practise in semi-analytic models of assigning concentrations 
to haloes according to just their mass and to some proba- 
bilistic acco unting of the dispersion shown in Fig. [7] As em- 
phasised bv lGao. Springel, fc White] l|2005l ). properties such 
as the clustering of haloes depend on the assembly history, so 
a full description of the correlation between formation his- 
tory and concentration may affect significantly the model 
predictions for the size and internal structure of a galaxy. 

A couple of prescriptions designed to predict halo con- 
centrations at z = from their mass accretion histories have 
been proposed recently (W02, Z03), and we use our simu- 
lation merger trees in order to compare them. We focus on 
the relaxed halo sample, since these have well-defined con- 
centration parameters. We also compare these methods with 
the simple prescription originally proposed by NFW. 
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Figure 12. The correlation between the halo characteristic den- 
sity, <5 C , and different definitions of formation time for relaxed 
haloes. The numbers along the top of each panel indicate the 
number of haloes in each bin. In the first three panels the forma- 
tion time is defined by reference to the most massive progenitor, 
and is set to be the redshift at which its mass was 1/2, 1/4 or 
1/10 of the final halo mass, Mo (see labels in each panel). In the 
final panel (labelled NFW) the formation time is defined as the 
redshift when half of the final halo mass is in progenitors more 
massive than 1/10 of Mo- The "natural scaling" 8 C <x (1 + z/) 3 is 
indicated by the solid line in each panel. 



3.4-2 Zhao et al. prescription 

Z03 differentiate two distinct phases in the MAH; one of 
early, fast accretion, followed by a slow-accretion period 
that lasts until the present. The transition between the two 
phases occurs at a characteristic redshift, z tp . This "turning- 
point" redshift may be used to estimate the concentration, 
assuming that the inner properties of the halo, such as the 
scale radius, r s , and its enclosed mass, are set at Ztp and 
vary weakly thereafter. 

This procedure is in principle straightforward to imple- 
ment in our simulations but we note that there are a sub- 
stantial number of haloes for which the distinction between 
the two accretion phases is not well-defined. In some cases, 
more than one phase of fast accretion seems to be present; 
in others, there is a single phase with no obvious turning 
point. This leads to ambiguities in the definition of zt p and 
its associated concentration that are not easily resolved and 
that affect a significant fraction of haloes. A similar worry 
applies to the Wechsler et al prescription, since eq. [7] is a 
poor approximation to the MAH of a significant number of 
systems. 



3.4-3 NFW prescription 

Finally, we consider NFW's proposal to identify the forma- 
tion redshift with the epoch when 50% of the halo is con- 
tained in progenitors more massive than certain fraction, /, 
of the final halo mass. NFW propose / = 0.01 in their orig- 
inal work in order to match the mass-concentration relation 
using the extended Press-Schechter formalism, but this frac- 
tion is too low to allow for an accurate estimate in N-body 
simulations. As a compromise, we adopt / = 0.1 for the 
results shown here. 



3-4-1 Wechsler et al. prescription 

W02 showed that the Mass Accretion History (MAH) of a 
halo's most massive progenitor, of mass Mo at redshift z — 0, 
may be approximated by a simple function, 

log 10 M (z) = log 10 Mo - az (7) 

i.e., by a straig ht line with slope —a in the plane log 10 M(z) 
vs. 2 (see also Ivan den Boschi (|2002l )). These authors then 
relate the parameter a to a formation time via at = 1/(1 + 
zi) — a/21n(10). Their Fig. 6 shows that this definition of 
formation time correlates well with the halo concentrations 
measured by B01 in their simulations and can be used to 
predict c at z = 0. 

Fitting this correlation, they find 

cw = Co ja u (8) 

with Co =4.1, the typical concentration of haloes forming at 
the present time. The implementation of this prescription in 
our simulations is straightforward, and we show the predic- 
tions in Fig. 1131 after recalibrating eq. [5] with cq = 2.26 in 
order to take into account our different definition of virial 
radius. 



3-4-4 Comparison between prescriptions 

Note that the three prescriptions described above are based 
on different features of the halo merger trees. While NFW 
looks at the mass spectrum of clumps containing half of the 
final halo mass, the other methods consider just the MAH of 
the most massive progenitor. The Z03 prescription depends 
on the slope of the scaling relation log 10 M s vs. log 10 r s in 
the slow accretion phase while the W02 recipe fits the whole 
MAH with a single slope. 

In spite of these differences, Fig. [13] shows that all 
three procedures yield concentrations that correlate rea- 
sonably well with the measured values. The rms scatter 
between prediction and measurement is indicated in each 
panel. It is smallest (marginally) for the NFW prescription, 
but even in this case it only reduces the scatter in the mass- 
concentration relation (Fig. [6} from o"i g M = 0.092 to ~ 0.077. 
Thus it only accounts for about 30% (|(of gM — o"NFw)/ <7 igMl) 
of the variance in the mass-concentration relation. 

The W02 prescription does similarly well by this mea- 
sure, but the slope of the c pre d-Cmeasured relation is a bit too 
shallow. The Z03 prediction has more scatter, but this is 
entirely due to a tail of haloes for which it predicts very low 
concentrations. We conclude that all three methods predict 
concentrations that correlate well with the measured values, 
but none of them is able to fully account for the scatter in 
the mass-concentration relation. 
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Figure 13. A compar ison of the measure d conc entrations with 
those predicted by the IZhao et alj l l2003ah (top) , IWechsIer et al.l 
j2002h (bottom left) and NFW (bottom right) prescriptions. The 
plotted contours enclose 65%, 95% and 99% of the haloes with the 
remaining 1% plotted as points. The cr-values indicated in each 
panel give the rms scatter in the prediction ((log 1 Q(c/c pret j)) 2 ) 1 / 2 , 
while <t igM is the corresponding rms scatter about the mass- 
concentration relation for the same set of haloes. 

4 SUMMARY 

We use the Millennium Simulation to examine the structural 
parameters of dark matter haloes formed in the ACDM cos- 
mogony. The large volume probed by the MS, together with 
its unprecedented numerical resolution, allow us to probe 
confidently the mass profiles of haloes spanning more than 
three decades in mass. Our main conclusions may be sum- 
marised as follows. 

• As in earlier studies, we find that the mass profile of dy- 
namically relaxed haloes are well approximated by the two- 
parameter NFW profile. We find that at least 1000 particles 
are needed in order to obtain unbiased estimates of the dis- 
tribution of halo concentrations and illustrate a number of 
potential pitfalls that arise from analysing poorly-resolved 
haloes, or from including in the sample haloes manifestly 
out of equilibrium. 

• We study the correlation between the NFW fit param- 
eters, which we express in terms of the halo mass and a con- 
centration parameter. These results extend previous studies 
to much larger halo masses than hitherto repo rted in the lit- 
eratur e. Combining our results with those of accio et al.l 
|2007l ). we find that a single power law reproduces the mass- 
concentration relation for over six decades in mass. These 
results are in reason able, albeit not perfect, agreement with 
the predictions of the lEke. Navarro, fc Steinmetz j200ll ) and 
NFW models. The model of iBullock et all l|200ll ) fails at 
large masses, and predicts concentrations at least a factor 
of ~ 2 too low for M ~ 10 15 h' 1 M© haloes. 

• The dependence of concentration on mass, while well 
established, is weak, and of equal importance is the broad 
scatter in concentration at fixed mass. The distribution of 



concentrations at given mass is well fitted by a lognormal 
function where both the mean and the dispersion decrease 
with increasing halo mass. These results allow us to estimate 
in detail the abundance of haloes with unusually low or un- 
usually high concentration, providing a well-defined predic- 
tion that may be used to interpret observations of objects 
of unusual density, such as highly-effective cluster lenses or 
galaxies with haloes of anomalously low density. 

• We find that, once unrelaxed haloes are excluded, there 
is no significant correlation betw een halo spin and concen- 
tration, contrary to the results of iBailin et al.1 (|2005h . 

• We have searched for several ways to account for the 
large dispersion in concentrations at given mass. The scat- 
ter in concentrations seems to arise largely due to variations 
in the formation time. We examined various plausible defi- 
nitions of the formation time, and find that concentrations 
are best predicted by formation times defined taking into 
account the collapse history of the full spectrum of progen- 
itors rather than the evolution of the single most massive 
progenitor. 

• We comp a re th e schemes of IZhao et al.l (|2003al ) . 
IWechsIer et ail (|2002l ). and a variant of the NFW prescrip- 
tion, and find that, while all three show a strong correlation 
between the predicted and measured concentrations, con- 
siderable scatter remains. In fact, none of these models is 
able to account for more than 30% of the intrinsic variance 
in the mass-concentration relation. It appears as if a large 
fraction of the scatter is truly stochastic or, else, dependent 
on aspects of the halo merger history that are not probed 
by these simple schemes. 
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